function [output] = fn_s_star(delta,R1,~,lambda,s_min,s_max,dist_d,d_min,d_max,phi)

% myfun = @(s) (R1-fn_rho1(s,phi))./(1-1./fn_rho2(s))-(1-lambda)*R1*zeta;
% 
% x0_min = 1.00001;
% x0_max = 1000*s_max;
% x0     = [x0_min, x0_max];
% 
% if myfun(x0_min)>0 && myfun(x0_max)>0
%     output = s_max;
% elseif myfun(x0_min)<0 && myfun(x0_max)<0
%     output = s_min;
% else
%     options      = optimset('Display','off');
%     [sol,~,flag] = fzero(@(s) myfun(s),x0,options);
%     output       = min(sol,s_max);
%     if flag ~= 1; disp('Error solving s_star'); end
% end

zeta  = fn_insured_deposits(delta,R1,dist_d,d_min,d_max);

C(1)  = phi;
C(2)  = -(R1*(1-(1-lambda)*zeta)-1+phi);
C(3)  = -(1-lambda)*R1*zeta;
sol   = roots(C);
s_sol = max(sol);

output = s_sol*(s_sol<s_max && s_sol>s_min) + s_min*(s_sol<=s_min) + s_max*(s_sol>=s_max);

end